rm(list=ls())
setwd("/Users/yuehou/Dropbox/KingReplication/Ongoing Work")
 

#############################
#############################
#### Loop with Beck-Katz ####
#############################
#############################

# NOTE: You must input the models (below) manually 

# This is a hack of zelig functionality, that allows for correct
# standard errors and coefficients, but also allows the use of
# non-zelig models (plm), etc
# It also allows the user to modify the imputed dataset by adding lags, etc


load("imputations23.Rdata")
m <- 8


acf1<-matrix(NA,m,6)

# Loads initial imputations into dataframe
for (j in 1:m){
	assign(paste("data_",j,sep=""), amelia.out9$imputation[[j]])
}


for (j in 1:m){
	data_temp <- get((paste("data_",j,sep="")))
	

	# Reframe as a panel dataset
	library(plm)
	p_newset <- pdata.frame(data_temp)
	
	p_newset2<-subset(data_temp,dhasminwage==1)  #just for minwage
	p_newset2<-pdata.frame(p_newset2)

	p_newset$log_oecd_gdp <- log(p_newset$oecd_gdp)
	p_newset$lag_oecd_gdp <- lag(p_newset$oecd_gdp)

# Create lag terms here
	p_newset$lag_ia5010 <- lag(p_newset$final_ia5010)
	p_newset$lag_govem_yh2<- lag(p_newset$govem_yh2)

	p_newset$lag_leftc<- lag(p_newset$leftc)
	p_newset$lag_leftc_2<- lag(p_newset$leftc,2)
	p_newset$lag_leftc_3<- lag(p_newset$leftc,3)
	p_newset$lag_leftc_4<- lag(p_newset$leftc,4)
	p_newset$lag_leftc_5<- lag(p_newset$leftc,5)

	p_newset$lag_leftgs   <- lag(p_newset$leftgs)
	p_newset$lag_leftgs_2 <- lag(p_newset$leftgs,2)
	p_newset$lag_leftgs_3 <- lag(p_newset$leftgs,3)
	p_newset$lag_leftgs_4 <- lag(p_newset$leftgs,4)
	p_newset$lag_leftgs_5 <- lag(p_newset$leftgs,5)

	p_newset$lag_arm_generos = lag(p_newset$arm_generos)
	p_newset$lag_minwag_jf = lag(p_newset$minwag_jf)
	
	

# Transform generosity measure

	temp1 <- p_newset$soc_exp1 <- (p_newset$arm_generos * p_newset$oecd_gdp) - (p_newset$lag_arm_generos*p_newset$lag_oecd_gdp)
	
	temp2 <- p_newset$lag_arm_generos * p_newset$lag_oecd_gdp
	
	temp3 <- p_newset$soc_exp <- 100 * (temp1 / temp2)
	
	p_newset$lag_soc_exp<- lag(p_newset$soc_exp) #we don't want this lag

	#assign(paste("as_",j,sep=""), p_newset$arm_generos)
	#assign(paste("as2_",j,sep=""), p_newset$lag_arm_generos)
	#assign(paste("as3_",j,sep=""), p_newset$soc_exp )


	# Oil covariate
	p_newset$oilshock <- 0
	p_newset$oilshock <- (factor(p_newset$year)==1973 | factor(p_newset$year) == 1972 | factor(p_newset$year) == 1971 | 	factor(p_newset$year) == 1970)



	
 ###############################################
   

# import model
   
    

    results<- plm(model="within",data=p_newset,soc_exp ~  lag_soc_exp+ vk_corp + vk_corp:leftgs + leftgs + vk_corp:lag_leftgs + lag_leftgs  + lag_leftgs_2 + lag_leftgs_2:vk_corp+ oecd_debt + arm_unemp + arm_open + arm_finop + oecd_gdpgr + oilshock) #2 lag
    
    
  
   results2 <- coeftest(results, vcov=function(x) vcovBK(x, cluster=c("group","time")))
 


   hkvalue<-seq(1,5,1)
   
   beta_interaction<- results$coeff["leftgs"]+hkvalue*results$coeff["vk_corp:leftgs"]
   
   
   covariance1 <-vcovBK(results,cluster=c("group","time"))["leftgs","vk_corp:leftgs"]
   
   
   se_interaction <- sqrt(results2["leftgs",2]^2+(hkvalue^2)*results2["vk_corp:leftgs",2]^2+2*hkvalue*covariance1)


results <- coeftest(results, vcov=function(x) vcovBK(x, cluster=c("group","time")))



# Output vectors
	assign(paste("resultsc_",j,sep=""),results[,1])
	assign(paste("resultss_",j,sep=""),results[,2])
	
	assign(paste("resultsinteraction_",j,sep=""),beta_interaction)
	
	assign(paste("resultsinterSE_",j,sep=""),se_interaction)

	



}





# Vectors of Results 	

coef <- 0
coef <- 0
coef_interaction<-0
se <- 0
se2 <-0
se_interaction<-0
se_interaction_2<-0


for (j in 1:m){
	coef <- coef + get(paste("resultsc_",j,sep=""))
	coef_interaction<- coef_interaction + get(paste("resultsinteraction_",j,sep=""))
	
	se <- se + get(paste("resultss_",j,sep="")) 
	se2 <- se2 + (get(paste("resultss_",j,sep="")))^2
	
	se_interaction<- se_interaction + get(paste("resultsinterSE_",j,sep="")) 
	se_interaction_2<- se_interaction_2 + (get(paste("resultsinterSE_",j,sep="")))^2	
}



coef <- coef/m
coef_interaction<-coef_interaction/m
se <- se/m
se2 <- se2/m
se_interaction_2 <-se_interaction_2/m
se3 <- 0
se4 <- 0


for (j in 1:m){
	se3 <- se3 + (get(paste("resultsc_",j,sep="")) - coef)^2
	se4 <- se4 + (get(paste("resultsinteraction_",j,sep=""))-coef_interaction)^2

}


se3 <- se3/(m-1)
se4 <- se4/(m-1)
final_se <- sqrt(se2 + se3*(1+(1/m)))
t <- coef/final_se

final_se_interaction<-sqrt(se_interaction_2 + se4*(1+(1/m)))
t_interaction<-coef_interaction/final_se_interaction




# Output Results
cbind(coef,final_se,t)

# Marginal Effects 
cbind(coef_interaction,final_se_interaction,t_interaction)





##### PLOT#####

vkrange=c(1,2,3,4,5)


plot(vkrange, coef_interaction,"l",lwd=2,ylim=c(-2.5,2.5),xlab="Corporatism",ylab="Marg. Effect of Left Gov",main="Welfare Generosity")

lines(lowess(vkrange,coef_interaction + 1.96*final_se_interaction),col="red",lwd=1,lty=3)

lines(lowess(vkrange,coef_interaction - 1.96*final_se_interaction),col="red",lwd=1,lty=3)

abline(h=0,col="blue",lwd=.5)








